load c:\manoj\projects\ace\Julia_W_new w w_climate Julia_W fday julia_euvac julia_solar_flux;
load c:\manoj\projects\ace\OMNI_ELEC_new;
load  c:\manoj\projects\longp\final_coefficients a_n b_n ;
load c:\manoj\projects\ace\substorm;
%w=(w_climate)*24.366*1e-6; %This copies the Julia drift (measured)-Julia drift (modeled) to w
%the multiplication converts drift in ms^-1 to electric field in Vm^-1
ace_fday = floor(ace_all(:,1));
phase_delay = 17;%minutes


a=a_n;
b=b_n;
 n_julia_prediction = 0;

% load c:\manoj\projects\substorm\selected_substorm;
% fday_substorm = selected_substorm;
% fday_substorm = non_selected_substorm;

%load c:\manoj\projects\ace\rand_substorm;
Julia_fday = floor(fday(:,1));

% kk = rand([1,4193]);
% fday_substorm = sort(kk)*2.1780e+003;

st = 14/24;
en = 20/24;
i_counter=0;
radius_epoch = 1/24;

for i = 1:length(fday_substorm),
    
    
    if fday_substorm(i)-floor(fday_substorm(i)) > st & ...
            fday_substorm(i)-floor(fday_substorm(i)) < en,
        
       [fday_int,ii] = intersect(Julia_fday,floor(fday_substorm(i)));
       if ~isempty(fday_int),
           
               
           
           if fday_substorm(i)-radius_epoch >= fday(ii,Julia_W(ii).k(1)) &...
               fday_substorm(i)+radius_epoch <= fday(ii,Julia_W(ii).k(end)),
           i_counter = i_counter+1;
           
           L = fday(ii,:) >= (fday_substorm(i) + 1e-5 ) - radius_epoch & ...
               fday(ii,:) <= (fday_substorm(i) + 1e-5 ) + radius_epoch;
           %1e-5 is to give a small offset to avoid fday_substorm at
           %exact time as a JULIA recording (hope so)
           
           w(ii,:) = w(ii,:) - nanmean(w); % remove the daily forcing
           
           SEG_JULI(i_counter,:) = w(ii,L);
           
           
           
            L = ace_fday == Julia_W(ii).fday;
            if sum(L) > 0,
            
            n_julia_prediction =  n_julia_prediction + 1;   
            ace_time = ace_all(L,1)+(phase_delay)/(60*24);%to delay ace data 
            ace_data = ace_all(L,2);
            L = isnan(ace_data);
            if sum(L) < 10 & sum(L) ~= 0,
            ace_data = interp1(ace_time(~L),ace_data(~L),ace_time);
            end;
           
            pred_julia =filter(b,a,ace_data);
            %pred_julia=pred_julia-mean(pred_julia);
            
             L = ace_time >= (fday_substorm(i) + 1e-5 ) - radius_epoch & ...
              ace_time <= (fday_substorm(i) + 1e-5 ) + radius_epoch;
          
            SEG_PRED_JULIA(n_julia_prediction,:) = pred_julia(L);
            SEG_IEY(n_julia_prediction,:) = ace_data(L);
            
            end;
           
           end;
       end;
    end;
end;

% plot((-radius_epoch*1440) + 2.5 : 5 :  (radius_epoch*1440) - 2.5 ,nanmean(SEG_JULI)*24.366*1e-3);
% hold on;
plot((-radius_epoch*1440) + 2.5 : 5 :  (radius_epoch*1440) - 2.5 ,nanmean(SEG_PRED_JULIA),'b');

